/******************************************************************************

                                BASICS.C

  Title:        General initialization package for SPEDEN.
  Author:       Hanna Szoke
  Date:         1/23/92
  Function:     This package contains set-up functions for all SPEDEN programs,
		involving the reading and interpretation of simple .inp file,
		and the echoing of its contents to the log.

******************************************************************************/

#include  "util.h"
#include  "cellparams.h"
#include  "dims.h"
#include  "symmetry.h"

  /************************************************************
  declarations of EXTERN variables in cellparams.h        
  ************************************************************/

real  aAxis, bAxis, cAxis ;	/* dimensions (nm or other) of unit cell	*/
real	angles[3] ;		/* 3 angles (degrees) of axes of unit cell */
real	crys_vol ; 		/* volfac() * aAxis * bAxis * cAxis	*/
float	incube_to_vox ;		/* conversion factor to el/voxel */
int	crystal_symmetry ;

  /************************************************************
  declarations of most EXTERN variables in dims.h               
  ************************************************************/

int     Nx, Ny, Nz ;            /* solution space dimensions */
float	floatNx ;
float	floatNy ;
float	floatNz ;
int     Np ;                    /* product of Nx*Ny*Nz */
int     Nptotal ;               /* Np */
int	Npextended ;		/* Nptotal + Nhr (high resolution)	*/
int	Nh, Nk, Nl ;		/* reciprocal space dimensions		*/
real  grid_spacing ;          /* gridding resolution in nm or other unit */
real	dr ;			/* actual average resolution */
real  eta ;                   /* factor used in setting up exp factors */
real	delfac ;		/* PISQ * eta * dr * dr	*/
real	input_res;		/* resolution input variable	*/
real	limit_rat ;    		/* 1/(dr*dr)	*/

int     useSig ;        	/* if TRUE, use sigmas in input fobs files */
int	HighRes ;		/* if TRUE, high-resolution processing on  */
float   min_voxel ; 		/* min. voxel value for a crystal (for ccg) */
float   max_voxel ; 		/* max. voxel value for a crystal (for ccg) */

 /************************************************************
  declarations of EXTERN variables in util.h               
  ************************************************************/

char	input_filename[MAXSTRING] ;
char	*allinp ;
char	nextline[MAXSTRING] ;
char	goodline[MAXLINES] ;
char	id[MAXSTRING] ;
char	units[MAXSTRING] ;
int	Nlines ;
int	HDFIN ;
int	HDFOUT ;

static	double	Na_error, Nb_error, Nc_error ;	/*** NOTE 16/2/05 ***/
static	float	mindens, maxdens ; 
static	char	record_info[MAXSTRING] ;  

static	void	dump_to_log() ;
static	int	fetchAllInput() ;
static	void	readCellInput() ;
static	void	get_res_params() ;
static	int	interpret_sym_group() ;
static	void	intercheck() ;
static	int	fetchNpdims() ;
static  void	fetchDensities() ;
static	void	basic_echo() ;
static	void	fetch_record_info() ;
static  void	get_IO_file_types() ;

     /************************************************************
     Read basic input: this is called by ALL Speden programs     
     ************************************************************/

void readBasicInput()    
{
     void	prepare_F_and_Finv() ;
     void	setupNusq() ;
     real	volfac() ;
     void	start_record() ;

     /************************************************************
     Dump the .inp file as-is to the .log file for reference.
     Fetch all the input, place it in a buffer, allinp.
     ************************************************************/

     check_exist(input_filename) ;

     sprintf(message, "\tReading file %s\n", input_filename) ;
     printTwice(message) ;

     dump_to_log() ;
     Nlines = fetchAllInput() ;

     /************************************************************
     Read input, identifying data in terms of keywords;
     a, b, c and input_res are always read.
     As of 9/2005, angle input is NOT read.
     Optionally eta is also read.
     ************************************************************/

     readCellInput() ;

     get_res_params() ;

     /************************************************************
     Identify the symmetry group (default: P1). 
     Check for consistency.        
     ************************************************************/

     crystal_symmetry = interpret_sym_group() ;
     intercheck() ;
		
     /************************************************************
     Define dimensions. 
     ************************************************************/

     Np = fetchNpdims() ;

     floatNx = (float) Nx ;
     floatNy = (float) Ny ;
     floatNz = (float) Nz ;

     Nptotal = Np ;
     Npextended = Nptotal ;
     Nh = Nx ;
     Nk = 2 * Ny ;
     Nl = 2 * Nz ;

     dr = sqrt(((aAxis/Nx)*(aAxis/Nx) + (bAxis/Ny)*(bAxis/Ny) + 
		(cAxis/Nz)*(cAxis/Nz))/3) ;

     Na_error = (Nx - (aAxis/dr)) / (aAxis/dr) * 100. ;
     Nb_error = (Ny - (bAxis/dr)) / (bAxis/dr) * 100. ;
     Nc_error = (Nz - (cAxis/dr)) / (cAxis/dr) * 100. ;

     limit_rat = 1 / (dr*dr) ;
     delfac = PISQ * eta * dr * dr ;

     /************************************************************
     Determine I/O file types.  Eventually, this will be HDF by 
     default.
     ************************************************************/

     get_IO_file_types() ;

     /************************************************************
     Use sigmas? Average solvent and protein voxel values?
     Minimum. maximum densities for ccg?
     Note: we call prepare_F_and_Finv "just to be on the safe 
     side" because it's needed in connection with pdb files.
     ************************************************************/

     setupNusq() ;
     useSig = fetchTFStatus("USESIG", TRUE) ;

     crys_vol = volfac() * aAxis * bAxis * cAxis ;
     incube_to_vox = crys_vol / Nptotal ;

     fetchDensities() ;
     min_voxel = mindens * incube_to_vox ;
     max_voxel = maxdens * incube_to_vox ;

     HighRes = fetchTFStatus("HIGHRES", FALSE) ;
     prepare_F_and_Finv() ;

     fetch_record_info() ;
     start_record();

     basic_echo() ;
}

     /************************************************************
     Dump all the .inp file as-is to the .log file.    
     ************************************************************/

void	dump_to_log()
{
     FILE       *fp ;

     fp = fopen(input_filename, "r") ;
     fprintf(fp_log,"Dump of %s follows:\n", input_filename) ;

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
          fprintf(fp_log, "\t%s", nextline) ;
     }

     fprintf(fp_log,"\nEnd of %s dump.\n", input_filename) ;
     fflush(fp_log) ;
     fclose(fp);
}

int	fetchAllInput()
{
     FILE       *fp ;
     int	N ;
     int	k, n ;
     int	nid ;

               /***********************************************
               Count lines of information in input file.
               ***********************************************/

     N = filelength(input_filename) ;

     if (N > MAXLINES)
         Fatal("Not enough space for input!  Please check your input file!") ;

     allinp   = (char *) e_malloc(sizeof(char)*MAXSTRING*N, 
		"fetchAllInput") ;

     for (n = 0; n < MAXSTRING*N; n++)
	  *(allinp+n) = '\b';

     for (n = 0; n < MAXLINES; n++)
	  goodline[n] = FALSE ;

     fp = fopen(input_filename, "r");

     N = 0 ;

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
          id[0] = '\0' ;
          if (net_strlen(nextline) > 1) {
               sscanf(nextline, "%s", id) ;
	       nid = (int)strlen(id) ;

               /***********************************************
               Identify data in terms of upper-case keywords.
               ***********************************************/

	       for (k=0; k<nid; k++)
		    nextline[k] = toupper((int) id[k]) ;

               sprintf((allinp+N*MAXSTRING), "%s", nextline) ;
	       N++ ;
          }
     }
     fclose(fp);

     return (N) ;
}

void readCellInput()    
{
     int	j, k ;
     char	setting[MAXSTRING] ;

     /**************
     Defaults 
     **************/

     t1 = t2 = t3 = t4 = t5 = t6 = 0 ;
     aAxis = bAxis = cAxis = 0 ;
     angles[0] = angles[1] = angles[2] = RIGHT ;
     strcpy(units, "nm") ;
     strcpy(setting, "") ;

     /* Note use of &t1 to capture single precision input! */

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "CELL") == 0) {
               sscanf(nextline, "%*s %g %g %g %g %g %g", 
			   &t1, &t2, &t3, &t4, &t5, &t6) ;
               aAxis = t1 ;
               bAxis = t2 ;
               cAxis = t3 ;
               /* 9/05: Angle input (t4 - t6) is unused */
	       goodline[k] = TRUE ;
          }
          if (strcmp(id, "UNITS") == 0) {
               sscanf(nextline, "%*s %s", setting) ;
               for (j=0; j<(int)strlen(setting); j++)
                    setting[j] = toupper((int) setting[j]) ;
               goodline[k] = TRUE ;
          }
     }

     /*********************************************
     Check for flawed or missing information. 
     *********************************************/

     if (strlen(setting) > 0) {
          if ((strcmp(setting, "ANGSTROM") == 0) ||
              (strcmp(setting, "A") == 0))
             strcpy(units, "A ") ;
          else if (strcmp(setting, "NM") == 0) 
             strcpy(units, "nm") ;
          else if ((strcmp(setting, "MICRON") == 0) ||
                   (strcmp(setting, "MU") == 0))
             strcpy(units, "mu") ;
          else
               Fatal(
          "Unknown units requested; please use Angstrom, nm or micron.") ;
     }     
}     

int	interpret_sym_group()
{
     char       symmetry_group[MAXSTRING] ;
     int	fetchSymop(char *) ;		/* packaged in crysutil.c */

/***
     int	k ;
     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "SYMMETRY") == 0) {
               sscanf(nextline, "%*s %s", symmetry_group) ;
	       goodline[k] = TRUE ;
               break ;
          }
     }
     for (k=0; k<(int)strlen(symmetry_group); k++)
         symmetry_group[k] = toupper((int) symmetry_group[k]) ;
***/

     /*********************************************
     In Speden, in fact we enforce P1 and disregard 
     SYMMETRY input.  Retrieve the symmetry system 
     and rotation-translation matrices.
     This is where Ncs (= 1 for P1) is established.
     *********************************************/

     strcpy(symmetry_group, "P1") ;		

     return (fetchSymop(symmetry_group)) ;
}

void intercheck() 
{

     /************************************************************
     The following checks are taken from International Tables for
     Crystallography, Volume A (1983), p. 13.
     There are no restrictions on triclinic crystals.

     Note that the b-unique setting is assumed throughout.  Note 
     too the assumption of hexagonal axes for all trigonal crystals.
     ************************************************************/

     if (strcmp(crystal_system, "MONOCLINIC") == 0) {
	if ((angles[0] != RIGHT) || (angles[2] != RIGHT)) 
	   Fatal("SPEDEN requires monoclinic crystals with b axis unique.") ;
     }		 
     else if (strcmp(crystal_system, "ORTHORHOMBIC") == 0) {
	if ((angles[0] != RIGHT) || (angles[1] != RIGHT) || (angles[2] != RIGHT)) 
	   Fatal(
	     "All angles must be 90 degrees for orthorhombic crystals.") ;
     }		 
     else if (strcmp(crystal_system, "TETRAGONAL") == 0) {
	if ((angles[0] != RIGHT) || (angles[1] != RIGHT) || (angles[2] != RIGHT))
	   Fatal(
             "All angles must be 90 degrees for tetragonal crystals.") ;
	if (aAxis != bAxis) 
	   Fatal(
             "Cell dimensions a and b must be equal for tetragonal crystals.") ;
     }		 
     else if (strcmp(crystal_system, "TRIGONAL") == 0) {
	if ((angles[0] != RIGHT) || (angles[1] != RIGHT) || (angles[2] != 120.)) {
	   Warning(    
      "Angles must be 90, 90 and 120 degrees for trigonal crystals.") ;
	   Fatal("Speden does not handle rhombohedral axes.") ;
         }
	if (aAxis != bAxis) 
	   Fatal(
              "Cell dimensions a and b must be equal for trigonal crystals.") ;
     }		 
     else if (strcmp(crystal_system, "HEXAGONAL") == 0) {
	if ((angles[0] != RIGHT) || (angles[1] != RIGHT) || (angles[2] != 120.)) 
	   Fatal(
            "Angles must be 90, 90 and 120 degrees for hexagonal crystals.") ;
         
	if (aAxis != bAxis) 
	   Fatal(
             "Cell dimensions a and b must be equal for hexagonal crystals.") ;
     }		 
     else if (strcmp(crystal_system, "CUBIC") == 0) {
	if ((angles[0] != RIGHT) || (angles[1] != RIGHT) || (angles[2] != RIGHT)) 
	   Fatal("All angles must be 90 degrees for cubic crystals.") ;
         
	if ((aAxis != bAxis) || (aAxis != cAxis)) 
	   Fatal("All cell dimensions must be equal for cubic crystals.") ;
     }		 
     else if (strcmp(crystal_system, "TRICLINIC") != 0) {
	sprintf(message, "Unknown crystal system - %s", crystal_system) ;
        Fatal(message) ;
     }		 
}

#define	SC_RESFAC	0.6 	/* see paper V	*/

     /************************************************************
     The following were CHANGED from long-standing old values on 
     6/23/00.  User input of grid_spacing was also dropped.   
     NOTE: Default res factors and etas lead to a pdb B value of
     about 11-1/2.
     ************************************************************/

#define	NEW_SC_ETA	0.8 	/* used to be 1.0	*/

void	get_res_params()
{
     int	k ;

     /* initializations */

     input_res = UNSET ;
     eta = UNSET ;

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "ETA") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               eta = t1 ;
	       goodline[k] = TRUE ;
          }
          else if ((strcmp(id, "INPUT_RES") == 0) || 
                   (strcmp(id, "RESOLUTION") == 0)) {
               sscanf(nextline, "%*s %g", &t1) ;
               input_res = t1 ;
	       goodline[k] = TRUE ;
          }
     }

     /* put it all together ... 
     no more old-style input of grid_spacing;
     new style is obligatory input_res only  */

     if (input_res == UNSET) 
	  Fatal("Missing input resolution!") ;
     

	  if (eta == UNSET) 
               eta = NEW_SC_ETA ;

	  grid_spacing = input_res * SC_RESFAC ;

     if ((eta < 0.5) || (eta > 0.9)) {
	  sprintf(message, "Eta (%g) is out of range (0.5 - 0.9).", eta) ;
	  Fatal(message) ;
     }
}


void     get_IO_file_types() 
{
     int	k, n ;
     char 	setting[6], usetting[6] ;

     /* temporary defaults, Nov. 2, 2005 */

     HDFIN = FALSE ;
     HDFOUT = FALSE ;

     for (k = 0; k < Nlines; k++) {
          strcpy(setting, "     ") ;
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%127s %5s", id, setting) ;

          if ((strcmp(id, "HDFIN") == 0)|| (strcmp(id, "HDFOUT") == 0)) {
               goodline[k] = TRUE ;

               for (n=0; n<6; n++)
                    usetting[n] = toupper((int) setting[n]) ;

               if (strcmp(id, "HDFIN") == 0) {
                    if (strcmp(usetting,"FALSE") == 0) 
                        HDFIN = FALSE ;
                    else if (strcmp(usetting,"TRUE") == 0) 
                        HDFIN = TRUE ;
                    else
                        goodline[k] = FALSE ;
               }
               else {
                    if (strcmp(usetting,"FALSE") == 0) 
                        HDFOUT = FALSE ;
                    else if (strcmp(usetting,"TRUE") == 0) 
                        HDFOUT = TRUE ;
                    else
                        goodline[k] = FALSE ;
               }
          }
     }
}

void	fetchDensities()
{
     int	k ;

     mindens = MIN_DENSITY ;	 /* default: 0 */
     maxdens = MAX_DENSITY ;	 /* default: 1.e10 */

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "MIN_DENS") == 0) {
               sscanf(nextline, "%*s %g", &mindens) ;
               goodline[k] = TRUE ;
          }
          else if (strcmp(id, "MAX_DENS") == 0) {
               sscanf(nextline, "%*s %g", &maxdens) ;
               goodline[k] = TRUE ;
          }
     }
}

void fetch_coefs(real *C1, real *C2)    /* Read C1 and C2 */
{
     int	k ;

     *C1 = *C2 = 1. ;

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "C1") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               *C1 = t1 ;
          }
          else if (strcmp(id, "C2") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               *C2 = t1 ;
          }
     }

     sprintf(message, "Using coefficients (%g, %g).", *C1, *C2) ;
     printTwice(message) ;
}

void	fetch_record_info()
{
     void	setup_records(char *);
     int	k ;

     strcpy(record_info, "history") ;

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "RECORD") == 0) {
               sscanf(nextline, "%*s %s", record_info) ;
	       goodline[k] = TRUE ;
          }
     }

     setup_records(record_info) ;
}

int	fetchNpdims() 

     /**********************************************************************
     fetchNpdims is called by all Speden programs that deal with gridded 
     arrays of electrons/voxel and Fourier transforms.  It bumps dx, dy and 
     dz to the closest product of powers allowed by the FFTs.  For certain 
     symmetries, we must impose further conditions on the dimensions.  
     These are expressed in maxdenom[i], i = 0, 2 and reflect the need to 
     have integer indices when applying a translational symmetry operation 
     on a given (x, y, z).   For most space groups, maxdenom[i] = 1 or 2, 
     but for P41, e.g., maxdenom[2] = 4 and other possible values are 3 & 6.

     The actual effective resolution is determined by averaging 
     aAxis/Nx, bAxis/Ny and cAxis/Nz.         
     **********************************************************************/

{
     int	getFftDim(float, int) ;
     float	dx, dy, dz ;

     dx = aAxis / grid_spacing ;
     dy = bAxis / grid_spacing ;
     dz = cAxis / grid_spacing ;

     Nx = getFftDim(dx, maxdenom[0]) ;
     Ny = getFftDim(dy, maxdenom[1]) ;
     Nz = getFftDim(dz, maxdenom[2]) ;

     return (Nx * Ny * Nz) ;
}

void	basic_echo() 
{
     fprintf(fp_log, "\n") ;
     fprintf(fp_log, "Unit cell measures %6.2f by %6.2f by %6.2f %2s\n",
             aAxis, bAxis, cAxis, units) ;

     fprintf(fp_log, "Alpha = %6.2f, beta = %6.2f, gamma = %6.2f degrees\n",
             angles[0], angles[1], angles[2]) ;

     fprintf(fp_log, "Scale factor for converting el/A^3 to el/grid pt is %g\n",
		     crys_vol / Nptotal) ;
     fprintf(fp_log, "Symmetry is %s.\n", SGname) ;
     fprintf(fp_log, "Input resolution is %g %2s.\n", input_res, units) ;

     /************************************************************
     Display and solution dimensions.
     ************************************************************/

     fprintf(fp_log, "\n") ;
     fprintf(fp_log, "Gridding resolution is %g %2s, eta is %g\n", 
		      grid_spacing, units, eta) ;
     fprintf(fp_log, "Input unit cell partitions are %4.1f by %4.1f by %4.1f\n",
		      aAxis/grid_spacing, bAxis/grid_spacing, 
		      cAxis/grid_spacing) ;
     fprintf(fp_log, "Actual unit cell partitions are %d by %d by %d\n", 
		Nx, Ny, Nz) ;
     fprintf(fp_log, "Average resolution in %2s is dr = %g\n", units, dr) ;
     fprintf(fp_log, 
     	"Partition errors (%%) with respect to the average resolution\n") ;
     fprintf(fp_log, "  in a, b and c are %g, %g, and %g.\n", 
		Na_error, Nb_error, Nc_error) ;
     fprintf(fp_log, "Run summaries will be written to %s.\n", record_info) ;

}

#define	JMAX	110	   /* dimension of goodval[] */
#define	JTOP	1024	/* max dimension (Nx, Ny or Nz)	*/

int 	getFftDim(float x, int M)	
			/* Expand x to the closest number that is an even 
			product of powers of 2, 3 and 5  and divisible by M.
			This is called by fetchNpdims() in basics.c.	*/

{
int goodval[] =       {   2,   4,   6,   8,  10,  12,  14,  16,  
			 18,  20,  22,  24,  26,  28,  30,  32,  
			 36,  40,  42,  44,  48,  50,  52,  54,  
			 56,  60,  64,  66,  70,  72,  78,  80,  
			 84,  88,  90,  96,  98, 100, 104, 108, 
		        110, 112, 120, 126, 128, 130, 132, 140,
			144, 150, 154, 156, 160, 162, 168, 176,
			180, 182, 192, 196, 198, 200, 208, 210,
			216, 220, 224, 234, 240, 250, 252, 256, 
			260, 264, 270, 280, 286, 288, 294, 300, 
			308, 312, 320, 324, 330, 336, 350, 352,
			360, 364, 378, 384, 390, 392, 396, 400, 
			416, 420, 432, 440, 448, 450, 462, 468,
			480, 486, 490, 500, 512, 1024} ;
   int Mgoodval[JMAX] ;

   int	i = 0, j ;
   int	imax ;
   int	closest_i ;

   for (imax = 0, j = 0; j < JMAX; j++)
	if ((goodval[j] % M) == 0 ) 
             Mgoodval[imax++] = goodval[j] ;

   if (x > JTOP) {
	sprintf(message, " %g is too large a dimension - maximum is %d", 
		x, JTOP) ;
        Fatal(message) ;
   }
   else if (x <= Mgoodval[0]) 
	return (Mgoodval[0]);

   while (x > Mgoodval[i])
      i++ ;

   closest_i = ((Mgoodval[i]-x) < (x - Mgoodval[i-1])) ? i : (i-1) ;

   return (Mgoodval[closest_i]) ;
}
